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COMPUTER PROGRAM FOR SOLVING 
COMPRESSIBLE NONSIMILAR- BOUNDARY- LAYER EQUATIONS 
FOR LAMINAR, TRANSITIONAL, OR TURBULENT FLOWS 
OF A PERFECT GAS 

By Joseph M. Price and Julius E. Harris 
Langley Research Center 

SUMMARY 

A computer program is described which solves the compressible laminar, transi- 
tional, or turbulent boundary-layer equations for planar or axisymmetric flows by an 
implicit finite -difference procedure. The program was used to obtain the solutions pre- 
sented in NASA TR R-368. Turbulent flow is treated by the inclusion of either a two- 
layer eddy-viscosity model or a mixing-length formulation. The eddy conductivity is 
related to the eddy viscosity by the static turbulent Prandtl number which may be an arbi- 
trary function of the distance from the wall boundary. The transitional boundary layer is 
treated by introducing an intermittency function which modifies the fully turbulent model. 
The intermittency function describes the probability distribution of turbulent spots and 
ranges from zero for laminar flow to vinity for a fully turbulent flow. 

INTRODUCTION 

A number of finite-difference methods are currently available for computing the 
development of compressible turbulent boundary layers. (See, for example, refs. 1 to 8.) 
The numerical methods used to solve the governing equations in these references are 
generally different, in particular those given in references 7 and 8; however, the results 
are similar when common eddy-viscosity formulations are used. References 1 to 7 use 
implicit or Crank Nicolson type differences to reduce the governing differential equations 
to finite-difference form, whereas explicit -type differences are used in reference 8. A 
coupled solution technique is used in reference 7 (see also ref. 9) which requires no iter- 
ation procedure. However, the methods used in references 1 to 6 require iteration since 
the momentum and energy equations are uncoupled. Another difference, of course, among 
the various methods is in the formulation of the eddy-viscosity and turbulent Prandtl 
number formulations used to model the turbulent flux terms appearing in the mean flow 
equations. 



This report describes a computer program developed to solve the compressible 
nonsimilar-boundary-layer equations for laminar, transitional, or turbulent flows of a 
perfect gas. The program was used to obtain the results reported in references 7 and 10. 
A coupled, implicit finite -difference procedure similar to that used for laminar flows in 
references 9 and 11 is used to solve the momentum and energy equations without iteration. 
The program will solve problems for two-dimensional or axisymmetric flow geometry for 
the flow of a perfect gas. Currently, power-law and perfect-air viscosity (Sutherland's) 
relations are included in the code; however, any perfect gas can be treated by inserting 
the correct viscosity temperature relations. Transverse curvature effects are included 
with the option of being neglected if desired. The equations are solved so that nonsimilar 
terms may also be neglected if desired. 

Options are provided for either a two-layer eddy-viscosity model or an arbitrary 
mixing-length formulation. The turbulent Prandtl number may be either a constant or a 
function of distance normal to the wall boundary. The transitional region of the flow is 
modeled through an intermittency distribution which modifies the fully turbulent eddy- 
viscosity models. Transition location and extent must be specified either from experi- 
mental data or correlation relations. The laminar equations are recovered by setting the 
intermittency to zero. 


SYMBOLS 


Values are given in both SI and U.S. Customary Units. The measurements and cal- 
culations were made in U.S. Customary Units. 


A damping function, 26v/uj 

A* = A\ij/v 


Aln,Bln,Cln,Dln,>i 

J 

A2n>®2jj,C2j^,D2jj7l 

E2n,F2n,G2jj J 


coefficients in difference equation (43a) and defined by equa- 
tions (B3) to (B9) 

coefficients in difference equation (43b) and defined by equa- 
tions (BIO) to (B16) 


Cf 


skin-friction coefficient. 


w 


1pu2 


CmljCmi defined in equations (B45) and (B46), respectively 
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specific heat at constant pressure 


Cp 

defined in equations (B36) and (B37), respectively 
Ey defined in equation (B39) 

F velocity ratio, u/ue 

^ml defined in equation (B29) 

^m2 defined in equation (B32) 

I 

I 

; Fy defined in equation (B40) 

G,H typical variables in the boundary layer (see appendix A) 

I Hi,H 2 ,H 3 ,. . coefficients defined by equations (B17) to (B28) 

h heat -transfer coefficient 

f i index used in grid-point notation (see eq. (41)) 

j flow index; j = 0 in planar flow, j = 1 in axisymmetric flow 

k grid-point spacing parameter (see eq. (41)) 

thermal conductivity 

k>p eddy conductivity (see eq. (15)) 

kj constant in eddy-viscosity model (see eq. (6)) 

k 2 constant in eddy-viscosity model (see eq. (7)) 

k 3 see equation (8) 

k 4 constant in intermittency function (see eq. (10)) 

L reference length 
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Lmi,Lml defined in equations (B34) and (B35), respectively 
I defined in equation (30) 

T mixing length (see eq. (13)) 

M Mach number 

m grid-point index in X-direction (see fig. 2) 

N number of grid points at each x-station (see fig. 2) 

Npr Prandtl number, CpjUy^k^ 

Npr t static turbulent Prandtl number (see eqs. (16) and (17)) 

Nst Stanton number, hy^^Cppu^ 

n grid-point index in Y-direction (see fig. 2) 

p(l),p(2),p(3) defined in equations (48) 
p pressure 

q(1),q(2),q(3) defined in equations (48) 
q heat -transfer rate 

R,Z body axes system with origin at the stagnation point, where Z is positive 

downstream and R is positive radially outward (see fig. 1) 

Rg unit Reynolds number, Ug/fg 

*^e,x Reynolds number based on x, Ugx/i/g 

^e,xt i Reynolds number at transition, UgX^ jy/t'g 

Re^g* Reynolds number based on displacement thickness, Ug6*y^ 
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Re ^0 Reynolds number based on momentum thickness, UgS/i^g 

Rg gas constant (see eq. (19)) 

r radial body coordinate measured normal to Z-axis (see fig. 1) 

rQ body radius (see fig. 1 ) 

S Sutherland's viscosity constant, 110.3° K (198.6° R) 

T static temperature 

Ty defined in equation (B41) 

adiabatic wall temperature 

Tmi,Tjji 2 defined in equations (B30) and (B33), respectively 

t transverse curvature term (see eq. (23)) 

u velocity component in X- direction (fig. 1) 

law of wall coordinate, u/u.^ 

Uj friction velocity, 

V transformed normal -velocity component (see eq. (26)) 

^ml defined in equation (B31) 

V velocity component in Y-direction 

V velocity component, v + 

X,Y orthogonal boundary-layer coordinate system with origin at stagnation point, 

where X lies along the body surface and is positive downstream and Y 
is normal to the body surface and positive outward (see fig. 1 ) 

Xi,X 2 ,...,Xs functions of grid-point distribution (see eqs. (A4) to (A 8 )) 
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X boundary- layer coordinate along X-axis (see fig. 1) 

x^ f end of transition (see fig. 1) 

xj I beginning of transition (see fig. 1) 

Yj,Y 2,. . .,Y0 functions of grid-point distributions (see eqs. (A12) to (A17)) 
y boundary-layer coordinate along Y-axis (see fig. 1) 

y+ law of wall coordinate, yu^/v 

y^jj match point for two-layer eddy-viscosity model 

z axial body coordinate (see fig. 1) 

a defined in equation (30) 

^ defined in equation (30) 

r streamwise intermittency distribution (see eq. (38)) 

y ratio of specific heats 

y transverse intermittency distribution (see eq. (10)) 

A* defined in equation (48g) 

Ax, Ay grid-point spacing, physical plane 

Axj transition extent, xj f - x^. j 

A^,A7 j grid-point spacing, transformed plane (see fig. 2) 

6 boundary -layer thickness 

6* displacement thickness 

* r°° 

incompressible displacement thickness, \ (1 - F) dy 

‘^O 
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eddy viscosity, 

eddy-viscosity function (see eq. (4)) 
defined in equation (B36b) 
eddy-viscosity function (see eq. (5)) 

transformed normal boundary-layer coordinate (see fig. 2) 

static -temperature ratio, Ty'Tg 

momentum thickness (see fig. 1) 

shock-wave angle (see fig. 1) 

defined in equation (40) 

molecular viscosity 

kinematic viscosity 

average kinematic viscosity 

transformed streamwise boundary-layer coordinate (see fig. 2) 

defined in equation (39) 

density 

exponent in power -law viscosity relation (see eq. (21)) 
shear stress 

local surface angle (see fig. 1) 

V by 


vorticity Reynolds number. 


V „ , maximum local value of y 
''^max 

stream function 

Subscripts: 

e based on boundary- layer edge conditions 

i inner region of turbulent layer 

m mesh point in | -direction (see fig. 2) 

n mesh point in rj-direction (see fig. 2) 

o outer region of turbulent layer 

r reference quantity 

s shock 

t total condition 

w wall value 

°o free stream 

Superscript: 

j flow index; j = 0 in planar flow, j = 1 in axisymmetric flow 

A prime on a symbol denotes a fluctuating component. 

A bar over a symbol denotes the time average value. 

A coordinate used as a subscript denotes the partial differential with respect to the 
coordinate. (See eqs. (Al).) 
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Figure 1.- Coordinate system and notation. 




PROBLEM DESCRIPTION 


This section presents the governing equations for compressible laminar, transi- 
tional, or turbulent boundary-layer flows together with the required boundary conditions. 
The eddy viscosity, eddy conductivity, transition location and extent, and transitional-flow- 
structure models are presented and briefly discussed; however, the reader interested in a 
detailed discussion of these models is referred to references 7 and 10. 


Basic Partial Differential Equations 

Governing equations .- The partial differential equations describing laminar, transi- 
tional, or turbulent compressible boundary -layer flows over planar or axisymmetric 
geometries are as follows (see refs. 7 and 12); 

Continuity 


Momentum 


P 


dy) dx 




Energy 



( 1 ) 

( 2 ) 


(3) 


where the conventional overbar notation for time mean-average variables has been 
dropped for brevity. The eddy-viscosity parameters ? and e are defined, respec- 
tively, as follows: 


and 



(4) 

(5) 


The intermittency-distribution parameter T is discussed in a subsequent section. (See 
eq. (38).) 

Eddy viscosity .- Options are provided in the coded program for selecting either a 
two-layer eddy-viscosity model (KODVIS = 1) or a straight-forward mixing-length model 
(KODVIS = 2). 
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Two-layer model .- The equations describing the two-layer model are as follows 
(see ref. 7): 



9u 

9y 


(ju)o " M ^2^e®inc>' 


where 


and 


D = 1 - exp 




VII 

VII 

(6) 

3 

A 

(7) 


( 8 ) 

(9) 


( 10 ) 


The boundary-layer thickness 6 appearing in equation (10) is defined as the distance 
normal to the wall boundary where u/Ug= 0.995. The empirical constants k 2 > k 3 , 

and k 4 are assigned values of 0.4, 0.0168, 0.0, and 0.78, respectively. Note that for 
k 3 = -1.0, the kinematic -viscosity term is removed from equation (8) (XT5 = -1.0). The 
empirical constants kj and k 2 can easily be treated as functions of some correlation 
parameter, say 6^, in order to account properly for low Reynolds number effects in 
hypersonic flows. (See ref. 13 for discussion of low Reynolds number effects.) The loca- 
tion of the boundary separating the two layers y^^ is determined from the continuity of 
eddy viscosity; that is, where 



( 11 ) 


Mixing-length model .- A mixing-length formulation is provided (KODVIS = 2) for 
those interested in utilizing experimental mixing-length distributions. (See ref. 13, for 
example.) The eddy-viscosity distribution across the boundary layer can be written as 
follows: 


i. = £ 

MM 9y 


( 12 ) 
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where the mixing length I may be written as 



( 13 ) 


Currently, the simplest possible formulation is provided in the digital code for 
follows: 



as 




(14) 


However, it should be noted that any functional variation can be utilized in the program. 

Eddy conductivity and static turbulent Prandtl number . - The eddy conductivity 
defined as 


v*t' 

= "CpP gr[./9y 


(15) 


is modeled as a static turbulent Prandtl number Npj. ^ as follows: 


NPr,t = 




(16) 


where 


Npr 


_ u*v* / 8 T/ 8 yV 

;^\9u/9yj 


(17) 


Any desired functional relation for Np^^t = ^^ 5 ) utilized in the digital code; 

three options are available. These options are (1) a constant, say Npj. ^ = 0.95, (2) an 
arbitrary distribution Np^^t = f^^ supplied in tabular form, and (3) the Rotta distribu- 
tion (see ref. 14) as follows: 


Npr,t = 0.45 



(18) 


The system of equations is closed by the addition of the perfect-gas laws and a 
viscosity-temperature relation. The perfect-gas law is expressed as 


P = pRgT (19) 

Currently, the digital code is written to include the Sutherland viscosity-temperature 
relation for air (IGAS = 1) 
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as well as the power-law expression 



where o = 0.647 for helium (IGAS = 2). 

Transformed plane. - The system of governing equations is singular at x = 0. The 
Probstein-Elliott (ref. 15) and Levy- Lees (ref. 16) transformation is used to remove this 
singularity as well as to reduce the growth of the boundary layer as the solution proceeds 
downstream. This transformation can be written as follows: 

l(x) = pgUeMe^o^dx (22a) 

where the parameter t appearing in equation (22b) is the transverse curvature term, 
defined as 


t = l+— (23a) 

^0 

or, in terms of the y-coordinate, as 

t = 1 + cos (p (23b) 

^o 

The relation between derivatives in the physical (x,y) and transformed (4 ,tj) coordinate 
system is as follows: 



= PeUeMe^o (^J + 
i/I? \Pe) 




(24a) 


(24b) 


Two new parameters F and 0 are introduced and defined as 



( 25 ) 



as well as a transformed normal velocity 


V = 


2g 

Pe'^e^e^o^ 



+ 


pvr|jt^ 

~w~ 


(26) 


The governing equations in the transformed plane can then be expressed as follows: 
Continuity 


W 

dn 


+ 2||f +F = 0 
9? 


(27) 


Momentum 


2{F^ + vH - s(f2 - o) > 0 

^ 9| 9tJ dT]\ drj ) ^ ' 


(28) 


Energy 


where 


2 |F^ * V|2 - ® ft2j * - aft2ie(|lf = o 

^ 9| 9tj 977 V Npj. 9 tj y \97] / 


I = 


a- 


(pm)© 


CpTe 




« _ 2^ (^e 


(29) 


(30) 


By using the viscosity relations (eqs. (20) and (21)) and the equation of state (eq. (19)), the 
parameter I can be written as follows: 


I = ^je 



(Air only) (31a) 
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(Power -law viscosity) 


(31b) 


I = (0)^"1-O 
where S = S/Tg. 

The transverse-curvature term can be written in terms of the transformed vari- 
ables as 

PeUe Jo p y 

The physical coordinate normal to the wall is obtained from the inverse transformation; 
namely, 


f’ e 

Pe^ero^ *^0 j 

The positive sign is used in equations (32) and (33) for axisymmetric flow over bodies of 
revolution (SIGN = 1.0), and the negative sign is used for flow inside axisymmetric ducts 
(SIGN = -1.0). 

The boundary conditions in the transformed plane are as follows: 

Wall boundary 


(33) 



F(^,0) = 0 
V(^,0)=V^(|) \ 

0(1,0) = ©w(^)^ 


or 



Edge conditions 



(34a) 


(34b) 
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The boundary condition at the wall for the transformed V component can be related to 
the physical plane as (see ref. 7) 


V ^ 


(35) 


Transition location .- Many parameters influence the location of transition. These 
parameters are discussed in an extensive review presented in reference 17. (See also 
ref. 7.) 

It is not currently possible to predict with assurance the location of transition for a 
general geometry; however, for some particular classes of flow such as those caused by 
sharp cones, empirical correlations are available which can be used with confidence pro- 
viding it is realized that a probable range of transition locations is being predicted and 
not an exact fixed point. (See ref. 7.) In the present digital code either the transition 
location (SST) or the stability index (SMXTR) must be specified; however, any correlation 
relation may be directly incorporated into the program if desired. 

Transition extent .- The assumption of a universal intermittency distribution implies 
that the transition-zone lei^h (transition extent) can be expressed as a function of the 
transition Reynolds number, In reference 18 it is shown, for the transition 

data considered, that the data are represented on the average by the equation 

^e,AX|. = ^^^e,x^ 


where Rg ^ ^Xj. f - x^. ^ . The location of the end of transition, 
obtained directly from equation (36) as follows: 


x^^f can then be 


^t,f = 


+ 5R, 


-1 



(37) 


where Rg is the local unit Reynolds number, vlq/Vq. 

In the present digital code, due to the lack of general correlations for the extent of 
transition, this quantity ^x^ f - xj. can be specified in one of two ways: (1) from equa- 
tion (37) (KTCOD = 1), or (2) from the specification of obt^-ined from experi- 

mental data (TLNGTH, KTCOD = 2). It should be noted that the digital code can be modi- 
fied to include any desired correlation or equation in place of equation (37). 

Intermittency distribution .- The parameter F appearing in equations (4) and (5) 
represents the streamwise intermittency distribution which models the turbulent spot 
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distribution in the transitional region; the parameter r is a function of the x-coordinate 
only and is defined as follows (see ref. 18): 

r(l) = 1 - exp(-0.412l^) (38) 


where 




X - xt,i 

X 


and 


(39) 


A = (x)j,^ - 

4 4 


(40) 


It should be noted that F = 0 for laminar flow, F = 1 for fully turbulent flow, and 
F ranges from 0 to 1 through the transitional -flow region. Equations (1) to (3) reduce to 
the classical laminar boundary-layer equations when F is set to zero. (See ref. 19.) 


Numerical Solution of the Governing Equations 

The system of governing equations (eqs. (19) to (21) and (27) to 29)) is parabolic and, 
therefore, can be numerically integrated in a step-by-step procedure in the streamwise 
direction. In order to cast the equations into a form in which the step-by-step procedure 
can be efficiently used, the derivatives with respect to ^ and tj are replaced by finite- 
difference quotients. The method of linearization and solution used in the analysis closely 
parallels that of references 9 and 11. 

Finite-difference mesh model .- It has been shown for laminar boundary layers that 
equally spaced grid points can be used in the normal coordinate direction. (See refs. 9 
and 11.) However, for transitional and turbulent boundary layers the use of equally 
spaced grid points is not practical because the fine-mesh size required to obtain accurate 
results near the wall boundary is inefficient for the entire boundary layer. The grid- 
point spacing in the rj-direction used in the program is such that the ratio of any two suc- 
cessive strips is a constant; that is, the successive Atj^ - coordinates form a geometric 
progression. In constructing the difference quotients, the sketch of the grid-point distri- 
bution presented in figure 2 is useful for reference. The dependent variables F, 0, and 
V are assumed known at each of the N grid points along the m-1 and m stations, 
but are unknown at station m+1. The and A ^2 values, not specified to be equal, 

are obtained from the specified x values (xin_l,Xui,Xj„^l) and from equation (22a). The 
relationship between the A? 7 ^-coordinates for the chosen grid-point spacing is given by the 
following equation (see ref. 1): 
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(i = l,2,3,...,N) (41) 


AtJj = (k)^"^ Atjj 

where k is the ratio of any two successive steps (XIQ, At 7 j^ is the spacing between the 
second grid point and the wall (note that the first grid point is at the wall boundary), and 
N denotes the total number of grid points across the chosen tj strip. The total thick- 
ness of the T} strip can then be expressed as follows: 




1 - (k)N-l 
1 - k 


(kj^l) (42) 


The selection of the optimum k and N values for a specified jjj^-coordinate depends 
upon the particular problem under consideration. The main objective in the selection is 
to obtain the minimum number of grid points with which a convergent solution may be 
obtained and thereby minimize the computer -processing time for each test case. The 
laminar boundary layer presents no problem since a k value of unity is acceptable; 
however, for transitional or turbulent layers, the value of k will be a number slightly 
greater than unity, say between 1.02 and 1.04. If transitional or turbulent flow occurs in 
a given problem, the laminar portion of the boundary layer is calculated with the value of 
k used for the turbulent region; that is, for a given problem, k is invariant. 

Difference equations .- Three -point implicit difference relations (see appendix A) 
are used to reduce the transformed momentum and energy equations (eqs. (28) and (29), 
respectively) to finite-difference form. The difference quotients produce linear difference 
equations when substituted into the momentum and energy equations provided truncation 
terms of the order A^^^ and Atj^.j A% are neglected, (it should be noted that 

the truncation term for d^F/drjl^ is of the order The resulting differ- 

ence equations may be written as follows: 

^^n^m+l,n-l ®^n^m+l,n ^ln^m+l,n+l ^^n®m+l,n-l 

+ Eln0m+i,n + ^ln®m+l,n+l = ^^n 


■^nEm+l,n-l ®^n^m+l,n 

+ E2n0m+l,n + F2n0m+l,n+l = (^^b) 

The coefficients Aljj,Bln,. . .,01^ and A2n,B2n,. . .,G2n (see appendix B) are functions 
of known quantities at stations m and m-1. It is important to note that equations (43) 
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are coupled through the dependent variables F and 0; however, the dependent vari- 
able V does not appear explicitly as an unknown at station m+1. The variable V is 

ft T? 

uncoupled from the system because of the particular way that the nonlinear terms V -r— 

OTj 

and (see eqs. (28) and (29), respectively) are linearized. (See eq. (A23).) 

Solution of difference equations .- The system of difference equations (eqs. (43)) 
represents a set of exactly 2(N - 1) linear algebraic equations for 2(N - 1) unknowns. 
The proper boundary conditions to be used with the difference equations are specified in 
equations (34). The 2(N - 1) linear algebraic equations may be written in tridiagonal 
matrix form; consequently, an efficient algorithm (Gaussian elimination) is available for 
simultaneous solution. 

The simultaneous or coupled- solution technique is presented in appendix B of ref- 
erence 9; however, because of differences between the present work and that presented in 
reference 9, the solution technique is discussed here in some detail. 

Because of the special form of equations (43), the following relations exist (see 
ref. 20): 


^m+l,n-l ^m+l,n-l ^m+l,n-l^m+l,n ^m+l,n-l®m+l,n 

®m+l,n-l ~ Qm+l,n-l + Qm+l,n-l^m+l,n Qm+l,n-l®m+l,n (44b) 

Next, equations (44) are substituted into equations (43) to obtain the following relations: 

®^m+l,n^m+l,n "*■ ®^m+l,n®m+l,n ~ ^^m+l,n “ ^^m+l,n^m+l,n+l 

^^m+l,n®m+l,n+l (45a) 

®^m+l,n^m+l,n ^2m+l,n®m+l,n “ ^2m+l,n " ^2m+l,n^m+l,n+l 

^lm+l,n®m+l,n+l (45b) 


where 


Blm+l,n - Blm+l,n + ^lm+l,nPS+i,n-l + ^lm+l,nQm+l,n-l 
®^m+l,n = ^^m+l,n + -^^m+l,nPm+l,n-l ^^m+l,nQm+l,n-l 
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(46c) 


^^m+l,n = ^^m+l,n ~ ^^m+l,n^m+l,n-l " ®^m+l,nQm+l,n-l 
^^m+l,n “ ®^m+l,n ■^^m+l,n^m+l,n-l ^^m+l,n^m+l,n-l (46d) 

®2m+l,n = ®‘2m+l,n '^2m+l,n^m+l,n-l ®^m+l,nQm+l,n-l (46e) 

^^m+l,n ~ ^^m+l,n “ ^^m+l,n^m+l,n-l ” ^2m+l,nQm+l,n-l (46f) 


The unknown values of F and 0 at station m+l,n are obtained from equations (45) 
as follows: 

^m+l,n = ^m+l,n "*■ ^m+l,n^ni+l,n+l + ^m+l,n®m+l,n+l (47a) 

®m+l,n “ Qm+l,n * ^m+l,n^m+l,n+l ■*■ *^m+l,n®m+l,n+l (47b) 

where 

^m+l,n ~ (^^m+ljn^lm+ljn " Elm+l,nG2m+l,n) ^m+l,n (48a) 

^m+l,n ~ (®^m+l,n^2m+l,n “ ®^m+l,n^^m+l,ii) ^m+l,n (48b) 

^m+l,n = (®^m+l,n^2jjj^.i^n - E2jjj+i n^^m+l,q) ^m+l,n (48c) 

^m+l,n “ (48c$ 


^m+l,n “ (®2m+l,n^^m+l,n " ®^m+l,n^2jjj^.i 

^m+l,n “ (®^m+l,n^^m+l,n “ ®^m+l,n^^m+l,n) ^m+l,n 

and 

^ _ 1 

' (BC+i,nE2^+l,n - B2*,+i,nEl^+l^„) 


(48e) 

(48f) 

(48g) 
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Next, equations (44) are rewritten as follows (where n = n+1) 


^m+l,n ^m+l,n + ^rn+l,n^m+l,n+l ^rn+l,n®m+l,n+l (49a) 

®m+l,n “ ^m+l,n "*■ ^m+l,n^m+l,n+l Qm+l,n®m+l,n+l (49b) 


The "no-slip" boundary condition j = 0^ is applied at the wall boundary to obtain 

the values of 1 where i = 1,2,3; that is. 


piili.i = pSIi.i " = » 


(50) 


The thermal condition at the wall boundary may be specified in one of two ways: 

(1) specified wall -temperature distribution, or (2) specified heat-transfer distribution. 
For a specified wall-temperature distribution it can be seen directly from equation (49b) 
that 


Qmll,l = ®m+l,l 
QSli,i = Qi^li,i = o[ 


(51) 


The case in which a heat-transfer distribution is specified presents a more difficult 
problem; however, this class of flows is often of interest. (For example, consider adia- 
batic flows.) 

The heat transfer at the wall boundary can be written in the transformed plane as 
follows (see ref. 7); 


_ ~Mr% ( Pe^e^c^^e"o ) , f 


(52) 


Then, for a specified value of the gradient of © can be obtained directly as 

follows: 


= -a / cuL \ l ^PrV^^ 


Npr^ 


(53) 
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For the grid-point spacing used in the analysis (geometric progression, see eq. (41)), the 
gradient of 0 evaluated at the wall, by using a 3 -point relation, is as follows: 




Qm+1,1 + - Qm+1,3 

k(l + k)Ar}| 


(54) 


Equations (53) and (54) then yield the following expression for 0m+l,l* 

_k(l^H^^/a©\ (1 + k)2 ^ 1 

m + 1,1 " 1 - (1 + k )2 1 - (1 + k )2 1 - (1 + k )2 


where (|^j is evaluated from equation (53). Equations (43) are next written at 

\®’’/m+l,l 

the m+1,2 point to obtain two equations in terms of ^ and ©m+l^n where 

n = 1,2,3. ^Note that = 0.^ The quantity Fi„+1,3 is next eliminated from these 

two equations to obtain one equation in terms of Fin+ij2 ®m+l,n where n = 1,2,3. 

The quantity ©m+l 3 is next eliminated through use of equation (55) to obtain the 
relation 




(56) 


where 




[(C2)(G1) - (C1)(G2)]^^1 [(C2MF1) - (C1)(F2)]„^, ,[>(1 . , 


^m+1,2 


(57a) 


,(2) |(C1)(B2) - (C2)(Bl)]^^i 2 


QJSil,! = 


^m+1,2 


(57b) 


_( 3 ) [(C1)(E2) - (C2)(E1)]^^^ 3 [(C1)(F2) - (C2)(Fl)]^„^i 


(57c) 


and 

Am+1,2 = [[(C2)(D1) - (C1)(D2)] + [(C2)(F1) - (C1)(F2)] [l - (1 + 2 
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By comparing equations (49b) and (56) it is observed that 


<‘ = 1 . 2 . 3 ) ( 58 ) 

which completes the desired boundary condition for the case of a specified heat-transfer 
distribution along the wall boundary. The temperature at the wall is obtained directly 
from equation (55) once ©ni+1,2 ®m+l,3 known. 

The quantities P^+l n Qm+l,n where i = 1,2,3 (see eqs. (49)) must first 

be determined across the boundary layer at the m+1 station where n = 1,2,. . .,N. 

These quantities are calculated by the following procedure: 

(1) Perform the following steps at the first grid point away from the wall (n = 2) : 

(a) Calculate Al 2 ,Bl 2 ,. . .,Gl 2 from equations (B3) to (B9). 

(b) Calculate A22,B22,. . .,G22 from equations (BIO) to (B16). 

(c) By using the results from steps (a) and (b) and the boundary conditions 

(eqs. (50) and (51) or (57)), calculate BlJ, B2|, El|, E22, GI 2 , 
and G2| from equations (46). 

(d) By using the results from steps (a) to (c), calculate and 

where i = 1,2,3 from equations (48). 

(2) The procedure outlined in step (1) is now repeated at grid point with n = 3 by 

using the results obtained at n = 2. This procedure is repeated until the 
entire boundary layer is traversed (n = N) and all values of ei. n and 
^m+1 n determined where i = 1,2,3 and n = 2,3,4,. . .,N. 

(3) By knowing the values of Pm+l,n Qm+l,n where i = 1,2,3 and 

n = 2,3,4,. . .,N, the values of and where n = N-l,N-2,. . .,2 

are calculated from equations (47). It should be noted that n 

®m+l,N specified e<^e boundary conditions (eqs. (34b)). The wall- 
boundary values of F and 0 are obtained from equations (34a) or equa- 
tion (55) for the case of a specified wall-boundary heat-transfer distribution. 
Before the computations can proceed downstream, the transformed velocity 
Vm+i,n must be determined across the boundary layer where n = 2,3,. . .,N. 
This requires the solution of the continuity equation. (See eq. (27).) 

Solution of continuity equation . - The continuity equation (eq. (27)) is solved numeri- 
cally for the N - 1 unknown values of V at station m+1. Equation (27) is integrated 
once to yield the following relation for Vm+l,n: 


24 




( 59 ) 


I d?7 

m +1 

where ^ represents the boundary condition at the wall V^. (See eq. (35).) The 

integral appearing in equation (59) is numerically integrated across the ? 7 -strip to obtain 
the N - 1 values of V. In the present program the trapezoidal rule of integration is 
used. 

Initial profiles .- Initial profiles for starting the finite-difference scheme are 
required at two x-stations since three-point differences are utilized. The initial pro- 
files at the stagnation point or line for blunt bodies, or near x = 0 for sharp-tipped 
bodies, are obtained by numerically solving the similar boundary-layer equations. (See 
eqs. (B47) to (B49).) The equations are solved by a fourth-order Runge-Kutta scheme 
with a Newton iteration method to modify the initial estimates of the gradients of F and 
© evaluated at the wall boundary. The N - 1 values of F, 0, and V obtained at the 
equally spaced N - 1 grid points are numerically redistributed to N - 1 grid points 
whose spacing is determined from equations (41) and (42) if a variable spacing is required. 
(As noted previously, variable spacing is required if transitional or turbulent flow occurs.) 
The second initial profile located at station m is assumed identical to the one located at 
station m-1. Any errors that might be incurred because of this assumption are mini- 
mized by using an extremely small value of A?; that is, an initial step size in the physical 
plane on the order of Ax = 1 x 10"® is used. The solution at the unknown station m+1 
is then obtained by the finite-difference method. Extremely small, equally spaced A|- 
steps are used in the region of the initial profile. The step size is increased after errors 
due to the starting procedvire have approached zero, that is, after 10 to 15 steps in A^. 

Evaluation of wall derivatives .- The shear stress and heat transfer at the wall are 
directly proportional to the gradient of F and © evaluated at the wall, respectively. 

By using G to represent a general quantity, where Gj^+i^l is the value of G evalu- 
ated at the wall, the four-point difference scheme used to evaluate derivatives at the wall 
is given as 


\8n ) " '^7^111+1,1 + Y3 Gui+i^ 2 + YgGn^+ 1^3 + Y^oGin+1^4 

\ '^m+ 1,1 

where the coefficients Y 7 ,. . .,Yjq are defined by the following relations: 

(1 + k + k2)2|k(l + k) - 1 ] + (1 + k) 

(1 + k) (1 + k + k2)k® At/j 


(60) 


(61a) 
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_(l+k + k 2 ) 

k2A,i 


(61b) 


Y - (l + k + k2) 
^ (1 + k)k^ Atjj 


(61c) 


and 


YlO = 


(l + k + k2)k^ Atjj 


(61d) 


For the case of equally spaced grid points in the tj - direction (k = 1), equations (61) 
become 


Y - 11 

6 Atj 

(62a) 

Y - 18 

8 6 A77 

(62b) 

Y - 9 

9 " 6 Atj 

(62c) 

V 2 

= 6 A, 

(62d) 


and equation (60) reduces to the familiar four-point relation; that is, 

PROGRAM DESCRIPTION 
General Discussion 

The program, written in FORTRAN IV for the CDC 6000 series computers, consists 
of three main programs, D2390, D23901, and D2401. Program D2390 computes the initial 
similarity solution to equations (B47) to (B49) with the points equally spaced in the t}- 
direction. Program D23901 takes this solution and redistributes the values of 77 geomet 
trically in subroutine GEOM, then interpolates to provide the solution over the geometri- 
cally spaced points. Program D2401 takes these data as a starting profile, reads other 
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input, and computes initial conditions. Steps are then taken down the body to solve the 
momentum and energy equations in finite -difference form, and the continuity equation is 
numerically integrated. Various boundary-layer parameters such as boundary-layer 
thickness, displacement thickness, momentum thickness, and skin-friction and heat- 
transfer coefficients are then calculated and the output is printed. Program D2401 uses 
the following subroutines: TURBLNT calculates the eddy viscosity, its derivatives, and 
the intermittency distributions required for the solution of transitional and turbulent flows; 
VARENT reads the variable -entropy input in tabular form, computes drs/dz, and prints 
the input and the derivatives; TABLE reads the body-geometry input, nondimensionalizes 
it, and, if necessary, distributes the values according to specified steps and computes 
derivatives; SETUP determines from input where profiles and wall values are to be 
printed; function INTEGT integrates by using the trapezoidal rule; INUNIT converts data, 
if necessary, to the U.S. Customary System of Units for computation and back to the 
International System for printout; FTLUP performs a second-order interpolation to find 
intermediate values from a tabular array (see appendix C). 

Descriptions, Flow Charts, and Listings of the 
Main Programs and Subprograms 

Main program D2390 .- The main program D2390 performs the locally similar solu- 
tion to the continuity, momentum, and energy equations ((B47), (B48), and (B49), respec- 
tively), by using a fourth-order Runge-Kutta technique with Newton's iteration method. 

The flow diagram of main program D2390 is as follows: 



i 
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Main program D23901 .- Main program D23901 takes the initial profile data from the 
program D2390 solution and redistributes the equally spaced tj values geometrically 
according to input-distribution constant XK and then interpolates to redistribute the 
profile to correspond with the new tj values. The flow diagram of the main program 
D23901 is as follows: 


( D23901 ) 

Read NAMELIST input 
from TAPE 9 



I 


Read profile j 

\ 

CALL GEOM 


Redistribute equally spaced ; 
.^77 values to geometrically^ ' 
spaced tj values 


CALL FTLUP 

^ Redistribute profiles 
\with new tj values/ 



(~Stop ] 
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Subroutine GEOM .- Subroutine GEOM redistributes the equally spaced tj values 
to geometrically spaced values according to input -distribution constant XK. The flow 
diagram for subroutine GEOM is as follows: 


GEOM 


I 



i 


Return 
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The program listing for subroutine GEOM is as follows; 
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Main program D2401 .- The main program D2401 controls the finite -difference solu- 
tion of the boundary-layer equations. It reads the initial profile data (which may come 
from D2390 or D23901) and other input, computes initial conditions, solves the momentum 
and energy equations in finite -difference form, numerically integrates the continuity equa- 
tion, calculates the boundary- layer thickness, displacement thickness, momentum thick- 
ness, and skin-friction and heat -transfer coefficients, and prints the output. The flow 
diagram of the main program D2401 is as follows: 


( D2401 ] 

I 

Initialize program 
constants and variables 


1 

/ Read NAMELIST input 


i 



38 


I 

© 






I 


(D 
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CALL SETUP^ 

Set up 
print control 


Calculate program 
constants 


CALL TABLE 



Read tabular input 
nd nondimensionalize, 
if necessary 
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Subroutine TURBLNT .- Subroutine TURBLNT calculates the eddy viscosity, its 
derivatives, and the intermittency distributions required for the solution of transitional 
and turbulent flows. The flow diagram for subroutine TURBLNT is as follows: 
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Subroutine VARENT .- Subroutine VARENT reads the variable -entropy input in tabu- 
lar form, computes drg/dz, and then writes the input and the derivatives. The flow 
diagram for subroutine VARENT is as follows: 
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Subroutine TABLE . - Subroutine TABLE reads tabular input for body geometry, non- 
dimensionalizes the input if necessary, distributes the values according to specified steps, 
and computes the derivatives. The flow diagram for subroutine TABLE is as follows: 
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The program listing for subroutine TABLE is as follows: 
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Subroutine SETUP .- Subroutine SETUP determines from the input where profiles 
and wall values are to be printed. The flow diagram for subroutine SETUP is as follows: 



( SETUP 1 

1 _ ^ 



i 

Determine stations to be 
printed as given by 
input data 

( 


t 


RETURN 
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Function INTEGT .- Function subroutine INTEGT integrates by using the trapezoidal 
rule. The flow diagram for function INTEGT is as follows: 


INTEGT 


Integrate by using 
trapezoidal rule 


T 

^ Return ^ 
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The program listing for function INTEGT is as follows: 
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Subroutine INUNIT .- Subroutine INUNIT converts International System dimensional- 
input data to the U.S. Customary System of Units for calculations in the program. The 
subroutine then converts the data back to the International System before output. The 
flow diagram for subroutine INUNIT is as follows: 
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USAGE 


The programs are run on the Control Data 6000 Series computer under the SCOPE 
3.0 operating system. The CPU time required for running all three programs is approxi- 
mately 0,003 second per mesh point. 

Array Dimensions 

Program D2401 uses the variable-dimension capability of the preprocessor 
installed at the Langley Research Center to enable the user to use a minimum amount of 
storage for each case. If this capability is not available at the user’s installation, the 
dimension statements at the beginning of program D2401 should be modified by inserting 
the following numbers in place of their equivalent designations: 

JK maximum number of steps in t 7 -direction plus 10 

JL = 1, if case is considering only constant entropy 

maximum number of steps in X-direction, if case is considering variable 
entropy 

JM maximum number of wall -value stations to be printed 

JN maximum number of profiles to be printed 

FORTRAN statements setting these values should also be inserted immediately following 
the NAMELIST statements in program D2401. 

Intermediate Data Storage 

The output for the initial solution found in program D2390 is written on TAPE 9 as 
well as on the output file. Program D23901 will then read the data from TAPE 9, redis- 
tribute the points geometrically, and write the redistributed solution in place of the origi- 
nal distribution on TAPE 9 as well as on the output file. Program D2401 will then read 
the redistributed solution from TAPE 9. If redistributing the points is unnecessary, 
executing program D23901 may be eliminated and the solution from D2390 will be read 
directly by program D2401. Generally, TAPE 9 will be a disk file to be used only for a 
current run. However, if many cases are to be run with the same initial solution, a 
physical tape can be requested so that D2390 and D23901 need not be rerun for each case. 
When using TAPE 9 as either a disk file or a tape file, it is automatically rewound at the 
beginning of D23901 and again at the beginning of D2401. 

Input Description 

Input for all programs is standard CDC NAMELIST. Program D2390 reads input 
listed under $NAM1 and copies these input data as well as the output data onto TAPE 9. 
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Program D23901 then reads these data from TAPE 9 as input. No other input is required 
for D23901. Program D2401 requires the input from TAPE 9 (written by either D2390 or 
D23901) and the data found listed under $NAM2. Subroutine TABLE (in program D2401) 
requires the data listed under $NAM3. If the case being considered is using variable 
entropy, then subroutine VARENT (in program D2401) requires the data listed under 
$NAM4. 

Dimensional input and output may be in either the International System of Units 
(KODUNIT = 1) or the U.S. Customary System of Units (KODUNIT = 0). The following 
listing of input and output data gives the units in the International System, followed in 
parentheses by the units in the U.S. Customary System. Where no units are given, the 
data are nondimensional. 

The $NAM1 input data for program D2390 are given as follows: 

HPR rj, increment for which values will be printed and stored on TAPE 9 

(Default = 0.1) 

XEND (see fig. 2) (Default = 10.0) 

H Runge-Kutta integration increment (Default = 0.01) 


PR 

Npr (Default = 0.7) 

XKK 

S/Tg (Default = 0.0) 

BETA 

/3 (Default = 0.0) 

ALPHA 

a (Default = 0.0) 

XO 

7]q (Default = 0.0) 

YO(l) 

Vw (Default = 0.0) 

YO(2) 

Fw (Default =0.0) 

YO(3) 

( (Default = 0.0) 

YO(4) 

©w (Default = 0.0) 

YO(5) 

(Default = 0.0) 

vnJv/ 
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YO(6) 


(Default = 0.0) 



XK k (H k = 1, then program D23901 does not change data from D2390.) 

(Default = 1.0) 

IGAS = 1 for air (Sutherland's viscosity constant) 

= 2 for power -law viscosity relation (Default = 1) 

VISCON constant in power-law viscosity relation, newton-sec/m2 (lb-sec/ft2) 

VISPOW a, exponent in power -law viscosity relation 

KODUNIT =0 if all dimensional input and output are in the U.S. Customary System 
of Units 

= 1 if all dimensional input and output are in the International System 
of Units (Default = 0) 


The $NAM2 input data for program D2401 are given as follows: 


XMA 


PTl 

p^ newton/m2 (lb/ft2) 

TTl 

T^^,ok (OR) 

WAVE 

shock-wave angle at tip of sharp body or stagnation point of blunt body 


(Default = 90.0), radians (degrees) 

XYl 



\PoO^=0 

XY2 



\ «>/x=0 

XY3 

(^e)x=0 

G 

y (Default = 1.4) 


R Rg, gas constant (Default = 1716.0), m2/sec2-OK (ft2/sec2-°R) 
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su 


S (Default = 198.6), (OR) 


PR 

PRT 

IBODY 

J 

W 

FT 

KOBE 

KODWAL 

lENTRO 

CONE 

lENDl 

A 

DS 

KODVIS 

SST 


Npr (Default = 0.72) 

Npr^t (IJefault = 0.9) 

= 1 for stagnation-point flows 
= 2 otherwise 

j 

= 0 if transverse curvature is neglected 
= 1 if transverse curvature is included (Default = 0) 

= 1.0 for nonsimilar solution 
= 0.0 for similar solution (Default = 1.0) 

= 1 if both laminar and turbulent profile values are defined for diagnostic 
reasons after flow is fully turbulent 
= 0 otherwise (Default = 0) 

= 1 for specified temperature distribution 
= 2 for specified heat-transfer distribution (Default = 1) 

= 1 for constant entropy 
= 2 for variable entropy (Default = 1) 

cone semiapex angle (Default = 0.0), radians (degrees) 

number of steps in X-direction 

reference length (Default = 1.0), meters (feet) 

initial step length in X-direction (Default = 0.01), meters (feet) 

= 1 for two-layer eddy-viscosity model 
= 2 for mixing-length model (Default = 1) 

x-location at which transition occurs (Default = 1.0E08), meters (feet) 


80 



SMXTR critical vorticity Reynolds number (Default = 1.0E08) 

TLNGTH xt f/xt i (Default = 2.0) 

! CORP coefficient in equation (38) (Default = 0.412) 

[ 

. CONSTNT transition model (Default = 0.0) 

I XTl ki (See eq. (6).) (Default = 0.4) 

XT2 A+ (Default = 26.0) 

I XT3 k 2 (See eq. (7).) (Default = 0.0168) 

XT4 k 4 (See eq. (10).) (Default = 0.78) 

! 

' XT5 k 3 (See eq. (8).) (Default = 0.0) 

PROINC incremental x value for which profile printouts will be made 
= 0 if only certain specified profile printouts will be made 
(Default = 1.0), meters (feet) 

PRNTESTC incremental x value for which wall -value printouts will be made 

= 0 if only certain specified wall-value stations printouts are desired 
(Default = 0.1), meters (feet) 

I IPRO number of specified profile printouts desired (other than those determined by 

PROINC) (Default = 0) 

PROVAL array of IPRO specific x values for which profile printouts are desired 
meters (feet) 

I IPRNT number of specified wall -value printouts desired (other than those determined 

by PRNTINC) (Default = 0) 

PRNTVAL array of IPRNT specific x values for which printouts are desired, 

I meters (feet) 
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NAUXPRO =1 if auxiliary profile printouts are desired (see output description) 
^ 1 otherwise (Default = 0) 


BLNGTH =0 if using constant step size in X-direction 

length of body if using variable step size in X-direction (Default = 0.0), 
meters (feet) 

NPUTYPE = 1 for dimensional input 

= 2 for nondimensional input (Default = 1) 

KODPRT = 1 for constant Npj. 

= 2 for Rotta distribution 

= 3 for tabular Np^. ^ = f(y/6) (Default = 1) 

NUMBl number of values read into PRTAR and GLAR arrays if KODPRT = 3 

PRTAR turbulent Prandtl number array, used only if KODPRT = 3 (NUMBl values) 

GLAR y/6 array corresponding to PRTAR, used only if KODPRT = 3 (NUMBl 

values) 


KTCOD =1 if transition extent is calculated from equation (37) 

= 2 if transition extent is read in as TLNGTH (Default = 2) 


The $NAM3 input data for program D2401 are given as follows: 

NUMBER number of values read into $NAM3 tables (Maximum = 100) 

L order of interpolation to be used for $NAM3 table (Default = 1) 

PE pressure-distribution array (NUMBER values), newton/m2 (lb/ft2) 

Z axial -coordinate array (NUMBER values), meters (feet) 

RMI body radial -coordinate array (NUMBER values) (Default = 1.0), meters 

TW wall temperature-distribution array (NUMBER values), °K (°R) 


(feet) 
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QW wall heat -transfer -distribution array (NUMBER values), 

watts/m2 (Btu/ft^-sec) 

RVWALD mass flux at wall array, (NUMBER values) (Default = 0.0), 

newton-sec/m3 (lb-sec/ft3) 

S x-station array corresponding to above table inputs (NUMBER values), 

meters (feet) 

SS array of incremental values between adjacent x stations for computation 

(Maximum = 1000), meters (feet) 

The first three values for SS must equal DS for the starting procedure; 
that is, SS(1) = SS(2) = SS(3) = DS. 


The $NAM4 input data for program D2401 are given as follows: 

NUMBER number of values read into $NAM4 tables (Maximum = 100) 

RRS array of radial coordinates of shock wave (NUMBER values) 

ZZS array of axial coordinates of shock wave (NUMBER values) 

Output Description 

The output for programs D2390 and D23901 consists of printing and the intermediate 
data on TAPE 9 as discussed earlier in this section. In program D2390, the $NAM1 input 
data are printed, followed by the initial profile consisting of the following values: 

Initial profile 

ETA = 7j 

YO(l) = V 

YO(2) = F 

YO(4) = 0 

Y0(5)=|§ 



83 


In program D23901 the output is printed in a form identical to that in D2390 except that 
the profile is redistributed. This same output is repeated in program D2401 for conve- 
nience. Next, the $NAM3 input data are printed. If the particular case is considering 
variable entropy, this is followed by the $NAM4 input data. 

Next, the $NAM2 input data and the $NAM constants which consist of the following 
values are given: 


PIO 



TIO 

Ttoo 

Tr 


G 

y, ratio of specific heats 


REY 

Pco^co^ 

free-stream Reynolds number 

^ oo 


RTl 

kilogram /m 3 (slug/ft 3) 


PI 

p^, free-stream pressure, newton/m2 

(lb/ft2) 

T1 

T^, free-stream temperature, (°R) 


R1 

p^, free-stream density, kilogram/m3 

(slug/ft 3) 

U1 

u^, free-stream velocity, m/sec (ft/sec) 

AAl 

a^, free-stream speed of sound, m/sec 

(ft/sec) 

TREE 

Tj., reference temperature, °K (°R) 


VISREF 

Pj., reference viscosity, newton-sec/m2 

(lb-sec/ft2) 

RIO 

Pt, °0 



Next, the profile values are printed according to the specifications in the input. 
These consist of the following: 
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Laminar -profile values 


ETA 

Y/YE 

U/UE 

T/TE 

TT/TTE 

CROCCO 

PT/PTR 

M/ME 

FZ 

TZ 

VORTREY 

XLMll 

Additional 

YPLUS 

UPLUS 

UDEF 


n 

y/Ye 

F 


0 


TtAt.e 

Tt,e " 

^ j., total pressure ratio 
M/Mg, Mach number ratio 

,n 

(x)m+i n’ Reynolds number 

^^^m+l,n 

(P#i)g 




values for transitional and turbulent profiles 

yuT 

V 

u 

Uj 

Ue - n 

1 + r, effective viscosity parameter 


VISEFF 



Auxiliary-profile values 


V 

GRAD(U/UE) 

GRAD(T/TE) 

FCl 

DAMP 

EPl 

EP2 

EP 

MKDEL 


V (see eq. (27)) 

FZ 

TZ 




3u 


9y 


/m+l,n 


1 - exp(- 




k 



m+l,n 



im+l,n 




‘5/m+l,n 


(see eq. (6)) 
(see eq. (7)) 


Next, the wall -value stations are printed according to the specifications in the input. 
These consist of the following: 

X spatial coordinate, meters (feet) 

XI I 

RAD rg, body radius (see fig. 1) 

Z axial coordinate of body (see fig. 1) 

BETA pressure-gradient parameter (see eqs. (30)) 

TRFCT r, intermittency distribution (see eq. (38)) 
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) 


t 

1 

K 

f 

ti 

I 

I 


RVWALD dimensional mass flux at wall (see eq. (35)) 

Pe^6^ 

REDELT , Reynolds number based on local displacement thickness 

Me 

P0^0® 

RETHET ■ 

^e 

PeUeX 

REX , local Reynolds number 

Me 

PE pressure, newton/m2 (lb/ft2) 

pj,u| 

T 

TE 7 ^, e(%e temperature, °K (°R) 

^r 

RE — , edge density, kilogram/m3 (slug/ft3) 

Pr 

UE — , edge velocity, m/sec (ft/sec) 

Ur 


ME Mg, e<^e Mach number 

MUE — , newton-sec/m2 (lb-sec/ft2) 

Pr 

DPEDX “s— , pressure gradient 

ox 

9Tg 

DTEDX temperature gradient 

ox 

9Ue 

DUEDX velocity gradient 

DLTAST 6*, displacement thickness, meters (feet) 


THETA 9 , momentum thickness, meters (feet) 

D/T shape factor 

U 

TAUD T^, wall shear stress, newton/m^ (lb/ft2) 

CFE - — skin-friction coefficient based on ec^e condition 

|pe4 
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CFW 

QSD 

HD 

NSTE 

NSTW 

NUE 

NUW 

SWANG 

ZSHK 

RSHK 

ITRO 

TW/TT 

RFTRUE 

ROUSE 

DSMXO 

XD 

YE 

UTAU 


skin-friction coefficient based on wall density 

|pw4 

heat transfer, watt/m^ (Btu/ft^-sec) 

— , heat-transfer coefficient, watt/m2-OK (Btu/ft2-sec-°R) 

“ ^aw 

— " ■ , Stanton number based on etfee condition 

Cp(pu)e’ 

— — , Stanton number based on wall condition 
^pPw^e 

Nusselt number based on edge condition 
Nusselt number based on wall condition 
local shock-wave angle, degrees 
axial coordinate of shock wave 
local radius of shock wave 

number of iterations performed for variable entropy 

^w 

Tt,oo 

T - T 

recovery factor 

^t " Tg 
^max 



/Pr%L 

/ Pr 


5g, boundary- layer thickness, meters (feet) 

It 

Ut = J—, m/sec (ft/sec) 

V P 
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PTR 

reference total pressure 

YMP 

^m 

P20 

^t.e 

9 


pA 

OMEGA 



\ Pr ) 


Sample Cases 

Two sample cases are presented in order to illustrate the input and output quantities 
in relation to the test conditions of the particular case being considered. These cases 
include laminar flow over a blunt axisymmetric body and laminar, transitional, and turbu- 
lent flow over a flat plate. 

Case 1 .- An example of laminar flow over a blunt, axisymmetric body is given in 
reference 21. The body is a spherically blunted, 25° half-angle cone. The wind-tunnel 
test conditions are as follows: 

M^ = 7.95 

Pt oo = 6-31 X 10® N/m2 

Tx „ = 7.83 X 1q2 °K 


The wall-boundary pressure distribution used in the numerical solution was obtained by 
the technique presented in reference 22. (See ref. 7 for comparison of numerical results 
with experimental data.) This case requires 75000g storage on the CDC 6000 computers. 
The listing of the variable -dimension data for this particular case is as follows: 

JK = 110 
JL = 1 
JM = 120 
JN = 20 
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The listing of the input data for case 1 is given as follows 


$NAMl 
HDP = 0. 1 , 

XE-N0=10.0* 

H=0.0l • 

PR=0.72, 

XKK*0. 1A085, 

BETA=0.5f 

ALPHA=0*0, 

X0=0.3» 

Y0( n=0.0,0.0f 0.77,0.38t0.:>8t0.0» 
XK=U0, 

IGAS^lf 

V/TSCnN = 0.7l73»=-C8t 
VISP1W=C.647, 

KnniiNiT=Of 

% 


$NAM2 

XMA=7.95* 

PTlsl3l760,0t 

TTI=*1410,0, 

WAVF=90.0* 

G=U4, 

R=l7l6.0f 

SU=198.6, 

PR=0.72f 

PRT*3.9, 

IRnnv=l , 

J = li 

W = 0f 

PT = UO, 

»<nr)P=0» 

KHOWAL^lt 

CHNF:=75.0t 
TE^’01 = 103, 

A=l.O, 

DS=0.a005f 

KOOVIS^l, 

SST = 0.1E4^09, 
FMXTR=0. lE+09, 
TLNGTH=2.0, 
CnPP=0.4L2t 
C0NSTNT=^O,0t 
XTI*0, 4, 

XT2=26.0» 

XT3=0,Ol68t 

X^4*0.78, 

XT5-0,0t 
ppoiNr. = o*i, 

PR^jTTN^ = •005t 
ippn=i » 

PPOVAL ( 1 )=0.345 , 
TPRNT=1 t 

PRNTVAU 11=0. 345t 

NADXPPO=Of 

PLNGTH=C.0» 

NPiJTYP6=lf 

KOOPRT=2t 

KTC00=2* 

% 


HMAM3 

^UJMc^CR = 77f 

L = l, 

PF{ 11 = 1150.4,1149.43, 1146. 54, 1141. 72f 1134. 98, 1126. 32. 1 1 15. 76 f 1103.29, 
1088.94,1072.74, 1054.72,1034.89,1013.3,989.997,965.019,938.401, 
910.389,880.655, 849.844,817.414,783.869,749.276, 713.361,676.169, 

6 33.401, 599. 196, 559. 265 , 5 17. 655,492 . 393,475 . 551 ,456. 972, 439.0 14, 
421.517,404.537,388.006,371.935,356.33 5,341. 208, 326. 529,312.31, 
298.539, 285.218,272.305,259.829,247.761,236. 107,231.92,231.678, 
231.161,229. 78,228.924,227.986,226.962,225.869,224.718,223.511, 
222.252,220.957,219.634,218.288,216.896,211.374,207.75,205.714, 
202.654, 201.503, 209. 7 29,223. 925,235. 348,247. 107,255. 308,251.082, 
264.994,264.131,262.866,262. 141,262.049, 

7( 1 )=0.0,0. 0000185 3, 0.00006 226, 0.0001354,0.0002 381,0.0003708,0.0005337, 
0.0C07273, 0.000952,0.001208,0.001497,0.001819,0.002175,3.002565, 
0.002991,0.003455 ,0.003957,0. 004498,0.005082,0.005708,0.00638, 
0.007 1 01,0.0 0 7 8 7 3,0.0 0 8 6 9 9,0.0 0 95 84, 0. 0 1 0 5 3 , Oi 01 1 55 , 0. 01265 , 
0.01401,0.01456,0.01512,0.01568,0.01624,0.0168,0.01737,0.01794, 
0.01652,0.0191 ,0.01968,0.02026,0.02085,0.02144,0.02204,0.02263, 
0.02324,0.02384,0.02445,0.02505,0.02606,0.02813,0.02921 ,0.03031, 
0.03143,0.03259,0.03378,0.03501,0.03628,0.03759,0.03895,0.04036, 

0. 04183,0.04843,0.0536, 0.05718,0.06457,0.08496,0. 1067,0.1268, 

0. 1474,0. 1699,0.1895,0.21,0.2511,0.292,0.3347,0.3799,0.4203, 

PMI ( 11 = 0.0,0.001243 ,0.00227 7,0.003356,0.004443,0.005546,0.006648, 
0.007751 ,0.008856,0.009962, 0.01107,0.01218,0.01329,0.01439, 
0.0155,0.01661,0.01772,0.01883,0.01994,0.02105,0.02216,0.02327, 
0.02437, 0.02 548 ,0.02659,0.02 769, 0. 0288 , 0 . 0299 ,0 . 031 16 ,0. 03 164 , 
0.03211,0.03257,0.03301,0.03343,0.03385,0.03425,0.03464,0.03502, 
0.03539,0.03575,0.03609,0.03643,0.03675,0.03707,0.03737,0.03766, 
0. 03795, 0.03823,0.0387,0.03966,0.04016,0.04068,0.0412,0.04174, 
0.0423,0.04287,0.04346,0.04407,0.0447,0.04536,0.04605,0.04913, 
0.05154,0.05321 , 0 .05665 , 0. 06616 , 0 . 07631 , 0. 08568 , 0.09528, 0.105 8, 
C. 1149,0. 1245,0.1436,0.1627,0.1326, 0.2037,0.2225, 

TW( 1) = 77<'540.0, 

9VWAL0( 11=77*0.0, 

St 1 > =0. 0,0.001243,0.00 227 8,0.00 336,0.00445 7,0.005563,0.006676,0.007797, 
0.008924,0.01006,0.0112,0.01236,0.01352,0.0147,0.01588,0.01709, 
0.01831,0. 01954, 0.02079,0.02207,0.02336,0.02469,0.02604,0.02742, 
0.02883,0. 03 02 9, 0.03179, 0.03 335,0. 0352,0.03594,0.03666,0.0 3 738, 
0.0381,0.03881,0.03951,0.04021,0.04091,0.0416,0.04228,0.04297, 

0. 04365,0.04433,0.04501,0.04 568,0.04636, 0.04 703,0.047 71 ,0.04836, 

0. 04948,0. 0517 7,0.05295,0.05416,0. 05541 ,0.05668,0.058,0.05935, 
0.06075,0.0622,0.0637,0.06525,0.06688,0.07415,0.07986,0.08381, 
0.09196,0. 1145,0. 1385,0. 1606 , 0 . 1834, 0 .2 082 , 0. 2298 , 0.2524, 0 . 2977 , 
0.3429,0.39,0.44, 0.4345, 

SSt 1)=10*.C005,50*.C01,4?*.01, 
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Case 2 .- An example of planar flow is presented in reference 23. The wind-tunnel 
test conditions were as follows: 

M^ = 2.8 

Pj ^ = 9.997 X 105 N/m2 
Tf oo = 3.11 X 1 q2 Ok 

= 0 (Adiabatic flow) 

The location of transition was not reported in reference 23; therefore, for the present 
calculations it was assumed that transition occurred near the sharp leading edge of the 
flat-plate model. This case requires 120000g storage on the CDC 6000 computer. The 
listing of the variable-dimension data for this particular case is as follows: 

JK = 310 

JL = 1 

JM = 750 

JN = 40 


The listing of the input data for case 2 is given as follows: 


tNAWl 

HPR=0.4, 

XENn=120.0, 

HrO.Ol, 

PP=0.7, 

XKK=o.qu, 

B?TA*0.C, 
n PH4=3 .n5, 
xn=o.Ot 

YOC .0,0.0, 0.4r>2,2. 43, -3.04, 0.0, 
XK = 1.02 , 

IG 4 S* 1 , 

VI«;CQN=0. 7173F-OP, 

VISP0W=O.647, 

K001INTT = 0, 

$ 
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j;mam 2 

XMA=2.8f 

PTl»21024*0, 

TTl=550.0t 
WAVP^O.Ot 
XYl=1.0t 
XY2=1.0f 
XY3=2. 8* 

G=1.4t 

R*1716.e, 

SU=l98.6f 

PP=0*7, 

PRTsO.95, 

mOY=2t 

J=0, 

W=0 f 
FT^l.O, 

K00f: = 0» 

KOOWAL=2t 

IPNT‘^n=i, 

CONiE^^O.O, 

lENDl^lOOOt 

A=l.O, 

DS=O.Ol , 

KOnviS-lf 

$ST=0,lE+09f 

SMXTR=25C0*C, 

TLN‘GTH-2*0t 

CnPP=0.412» 

CONSTNT=0.0, 

XTl:=0*4f 

XT2=26*0t 

XT3=0.0l68f 

XT4*D,78t 

XT5=0.0, 

poniNC=2,OT 

PPNTINC=0*1» 

!PPO=14, , 

PRPVAL (I)=0.1f0*15f0.2f0«25t0.3»0*35»0*4»0.45»0.02f0*05t0.08f06t.07tl*t 
IPRNT=13, 

PRNTVAL(l)=2*Of4.0»C.2T0.25f0.3f0.35f0.4f0.45f0.02f0.05»0.08».06».07f 

NAiJXPPn = Ot 

8LNGTH=0.0t 

NPUTYPP=1, 

KonpPT =:2 ♦ 

KTCOO^l t 
$ 


SNAM3 

NUMREP = 41 » 
l = lt 

PE( n = 41* 774.6986 I 58056 3 f 

7U ) =0. Of C. 5 ♦ 1. Of 1.5 1 2. Of 2. 5 » 3. Of 3.5 f4. 0,4.5 f 5. Of 5. 5, 6. Of 6.5 1 7. Of 7*5 f 

8.jf8.5f9.0f9.5fl0.0fl0.5fll.0fll.5fl2.0fl2.5fl3.0fl3.5f 14.0,14.5,15.0, 

15. 5, 16. Of 16. 5 ,17. Of 17. 5 ,13. Of 1 8. 5, 19. Of 19. 5, 20. Of 

PMT( U = 4 l* 1 . 0 f 

TW( U =41*520.4439, 

ow( n=4i*o.Of 
RVWALK l)=4l*0.0f 

S(l) = 0. Of 0. 5,1. Ofl. 5,2. 0,2. 5 ,3. 0,3. 5 14. Of 4. 5, 5. 0,5. 5,6. 0,6. 5,7. 0,7. 5, 

8. 0. 8. 5. 9. 0.9. 5. 10. 0.1 0.5. 11. 0.11.5. 12. 0.12. 5. 13. 0.13. 5.14.0. 14. 5. 1 5.0, 

15.5. 16. 0. 16. 5. 17. 0.17. 5. 18. 0.1 8. 5. 19. 0.19. 5. 20.0, 

SS( 1) = 1000*.01, 

% 


Sample outputs for case 2 for regions where the flow is laminar and turbulent are given as follows 
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APPENDIX A 


DIFFERENCE RELATIONS 


Three-point implicit difference relations are used to reduce the transformed 
momentum and energy equations (eqs. (28) and (29)) to finite-difference form. It is 
assumed that all data are known at the solution stations m-1 and m. (See fig. 2.) 
Then, it is possible to obtain the unknown quantities at the grid points for the m+1 sta- 
tion. In the subsequent development the notations G and H are utilized to represent 
any typical variable. 

Taylor-series expansions are first written about the unknown grid point (m+l,n) in 
the 4 "direction as follows: 


^m,n ~ ^m+l,n ~ 


m+l,n 


^^2/ N 




i°iu) 


m+l,n 


(Ala) 


and 


Gm-l,n = Gm+l,n ‘ ^ 

6 \ ?^Vm+l,n 


(±ii^2Y 


(^^^)m+l,n 


(Alb) 


where subscript notation has been utilized to denote differentiation; that is, s — , 

and so forth. 

Equations (Ala) and (Alb) can be solved to yield 




HA 


XiG^^l n'^2^m,n + ^3^m-l,n ^ A^2 (^^1 + A^2) , 


m+l,n 


2A|2 




(A2) 


and 


r ^ r - r < ^^2 L , ^^2\ ^ 

^m+l,n - ^4^m,n ' ^5^m-l,n + 2 \ ^ A|i/ 


(A3) 


Terms of the order of 
errors of the order of 


A|j A^2> smaller, are neglected. This produces truncation 
A^l A |2 instead of A^2 reference 9 where two-point 
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APPENDIX A - Continued 


difference relations are used. The Xi,X 2 ,. . .,X^ coefficients appearing in equa- 
tions (A2) and (A3) are defined as follows: 

+2A^2 

1 “ A^l + A^2 


(A4) 


Xo = 2 


X3 = 2 


and 


X4 = 


X5 = 


+ A^2 

A^2 

+ AI 2 ) 

A^l + A|2 


A|2 

All 


(A5) 

(A6) 

(A7) 

(A8) 


Taylor-series expansions are next written about the unknown grid point (m+l,n) in 
the Tj-direction as follows: 






(A9a) 


and 




Ailn-l , 


(A9b) 


Equations (A9a) and (A9b) can be solved to yield 




(AlO) 


and 


(^] ~ ^4Gm+l,n+l “ ^5Gm+l,n ' ^6^m+l,n-l o ^lyrir) * 

v^/m+l,n 


(All) 
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APPENDIX A - Continued 


The Yj,Y 2 ,. . .,Y 0 coefficients appearing in equations (AlO) and (All) are defined as 
follows: 


Yi = 7 r 

A77^(A7Jn + 

Y 2 

2 A77„A77„_i 

Y = 2 

Y ^^n-1 

4 Arjn(AT7„ + Atj^.i) 

Y - ^^n-1 - ^^n 

^n„A77n-l 

and 

Y 


(A12) 

(A13) 

(A14) 

(A15) 

(A16) 

(A17) 


For the case of equally spaced grid points in the and t}- coordinates, equa- 

tions (A4) to (A 8 ) and (A12) to (A17) reduce to the following relations: 


Xi = 3^ 
Xg=4 
X3 = l > 
X4 = 2 

X5 = l> 


(A18a) 


and 
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APPENDIX A — Continued 



Y2 = 2Yi 


*5 




Y4 


1 

2Ar] 


Y5 = 0 


Y6 = Y4 

J 


(A18b) 


where and A 77 represent the spacing between the grid points in the and 
Tj -coordinates, respectively. 

Equations (A2), (A3), (AlO), and (All) can then be written for constant grid-point 
spacing as follows: 



'm+l,n 


^^m-fl,n ~ ^^m,n ^m-l,n 
2A4 



+ • . . 


(A19) 


Gm+l,n = 2Gm,n " Gm-l,n + + 


(A20) 


/82g\ ^m+l,n+l ~ 2Gm+l,n ^m+l,n-l 


A7?2 


+ • • • 


(A21) 


and 


/9G\ _ ^m+l,n+l “ ^m+l,n-l Atj2 ^ 

Wm.l,n= • • 


(A22) 


Equations (A19) to (A22) are recognized as the standard relations for eqvially spaced grid 
points. (See, for example, ref. 11.) 
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APPENDIX A - Concluded 


Quantities of the form that J^>pear in the governing equations must be linear- 

ized in order to obtain a system of linear difference equations. Quantities of this type 
are obtained from equations (A2) and (A3). 


The procedure used to linearize nonlinear products such as 
as that used by Flugge-Lotz and Blottner (ref. 9) and is as follows: 


W A^/ 


is the same 




(A23) 


where the terms 




and 




, , are evaluated from equation (All), but at the 

\?^/m,n 

known station m. By equatii^ G to H in equation (A23), the linearized form for 
quantities of the type obtained; that is, 


where 



(A24) 


The preceding relations for the difference quotients produce linear -difference equa- 
tions when substituted into the governing differential equations (eqs. (43)) for the conser- 
vation of momentum (eq. (28)) and energy (eq. (29)), respectively, since terms of the order 
of (A^)^ are neglected. 
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APPENDIX B 


COEFFICIENTS FOR DIFFERENCE EQUATIONS 

Equations (43) are the difference equations used to represent the partial differential 
equations for the conservation of momentum and energy, respectively. These equations 
are repeated for convenience as follows: 

^^n^m+l,n-l ®^n^m+l,n ^ln^m+l,n+l + ^^n®m+l,n-l 

+ EljiOm+i^n + ^^n®m+l,n+l “ (®1) 

•^2n^m+l,n-l ®2n^m+l,n ^^nFm+l^n+l 

+ E2nOni+l,n + ^n®m+l,n+l = ^2^ (B2) 

These equations are obtained from equations (28) and (29) and the difference quotients 
are presented in appendix A. The coefficients Aljj, Bljj, and so forth, in equations (Bl) 
and (B2) are functions of known quantities evaluated at stations m and m-1. (See 
fig. 2.) Therefore, equations (Bl) and (B2) can be solved simultaneously without iterative 
procedures. The coefficients Al^, Bl„, and so forth are as follows: 


Ain = Y 3 H 3 - 

(B3) 

Bl„ = XiH, - Y2H3 - Y5H11 + H5 

(B4) 

«n = Y 1 H 3 * Y4HH 

(B5) 

Mn = -Y6H4Fy 

(B6) 


(B7) 


(B8) 
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APPENDIX B - Continued 


Gin = HlFm2 + H4TyFy 

(B9) 

A2n = -2Y0H8 Fy 

(BIO) 

B 2 n = ^A 2 „ 

(Bll) 

a 

g 

1 

II 

csi 

0 

(B12) 

D2n = Y3H10 - Y0H12 

(B13) 

E2n = XiHi - Y 2 H 10 - Y5H12 

(B14) 

F2n = YiHio + Y4H12 

(B15) 


and 

02„ = HiT„2 + * H9(Ty/ 

The coefficients Yi,Y 2 ,. . .,Y 0 and Xj,. . .,X 5 are functions of the grid-point 
spacing and are defined in equations (A12) to (A17) and (A4) to (A 8 ), respectively. The 
coefficients . -,^12 defined as follows: 


Hi = ^m+l^ml^ 

(B17) 

H 2 = Yjni - Ltnl(^mlGml + ^ml^ml) 

(B18) 

H3 = -EmlBmlGml 

(B19) 


(B20) 

H5 “ ^m+l^ml 

(B21) 
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and 


^6 = -^m +1 

(B22) 


(B23) 


(B24) 


(B25) 


(B26) 

Hii=H 2 + H4Ty 

(B27) 

Hi 2 = H7 + 2 H 9 TY 

(B28) 


The undefined quantities appearing in equations (B17) to (B28) are defined as 
follows: 


^ml ~^4^m,n " ^5^m-l,n 

(B29) 

“ ^4®m,n “ ^5®m-l,n 

(B30) 

Vml=X4V^,n-X5Vm-l,n 

(B31) 

Fm2 = X2FJJJ n - XaFjjj.j 

(B32) 

Tm2 =X2®m,n ' X3©m-l,n 

(B33) 



Tml + (^) 

Wm+1 

(Air only) (B34a) 
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where 


and 
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a-1 


(Power law) 


(B34b) 


Lml = 


Lml 

fe) ' 


2Tml 

(t") 


LV^e/m+1 J 


(Air only) (B35a) 


^ml - ~ l)('^ml) 

(Power law) 

(B35b) 

Eml = 


(B36a) 

\ ^m-l,n ^m,n ^m+l,n 


(B36b) 

. fev)ni+l,n 
Eml- ^ 


(B37) 

Ey - ^4^m,n+l " ^5^m,n " ^6^m,n-l 

(See eq. (All)) 

(B38) 

Ey = Y4€m^n+1 ■ ^5^m,n " ^6^m,n-l 


(B39) 

Ey ~ ^4Em,n+l “ ^5Em,n ” ^6Em,n-l 


(B40) 

Ty = - YgOjji^n “ ^6®m,n-l 


(B41) 


(See eqs. (30)) 

(B42) 


“m+1 = 


fu^\ 

— I 

\Te/ 

^ ®^m+l 


(B43) 
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The transverse-curvature terms are contained in the quantities C„il and Cjjjj 
which appear explicitly in the H2, H3, H7, Hg, and Hg coefficients. The transverse- 
curvature term in the transformed plane (see ref. 7) may be written as follows: 


t2j = 1 , 2coj(W)^ cos 0 rn Q 

Jrt 


(B44) 


where t represents the ratio r/r^ and is a known quantity for the N-1 grid points at 
station m-1 and m. Then, the extrapolated values at m+l,n are obtained as follows 
where the parameter C is used to represent t2j: 


Cml=X4 C„„-X5 C„.i_„ (B45) 

^ml = Y4^m,n+1 " ^5^m,n “ Y6Cm,n-l (B46) 


Two quantities (symbols) as of now remain undefined. These are the code sym- 
bols FT and W which appear in equations (B17) and (B44), respectively. The code 
symbol W appearing in equation (B44) is used either to retain or neglect the transverse- 
curvature terms for axisymmetric flows; that is, W = 1 or 0, respectively. For planar 
flows, the transverse -curvature term does not appear since j equals 0. The code 
symbol FT (flow type) appearing in equation (B17) is used either to retain or neglect the 
nonsimilar terms in the governing differential equations; that is, FT = 1 or 0, respec- 
tively. If FT is assigned a value of unity, the solution to the nonsimilar equations 
(eqs. (27) to (29)) is obtained. If FT is assigned a value of zero, the locally similar 
solution is obtained; that is, the following system of equations is solved: 


Continuity 



F = 0 


(B47) 


Momentum 





(B48) 


(B49) 
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The governing equations for the locally similar system are obtained from equa- 
tions (27) to (29) by neglecting derivatives of the dependent variables F, 0, and V with 
respect to the streamwise coordinate The capability of obtaining locally similar solu- 
tions is desirable in that for a given test case the locally similar and complete nonsimilar 
solutions can be obtained for the identical program inputs and numerical procedures. 
Consequently, the effects of the nonsimilar terms on the boundary-layer characteristics 
for a particular case can be determined by a direct comparison of the results obtained for 
solutions for FT = 1 and 0, respectively. 



APPENDIX C 


LANGLEY LIBRARY SUBROUTINE FTLUP 


Language : FORTRAN 

Purpose : Computes y = F(x) from a table of values using first- or second-order interpolation. 

An option to give y a constant value for any x is also provided. 

Use : CALL FTLUP(X, Y, M, N, VARI, YARD) 

X The name of the independent variable x. 

Y The name of the dependent variable y = F(x). 

M The order of interpolation (an integer) 

M = 0 for y a constant. VARD(I) corresponds to VARI(I) for 
1=1,2,. . ., N. For M = 0 or Nil,y = F(VARI(1)) for any value of x. 

The program extrapolates. 

M » 1 or 2. First or second order if VARI is strictly increasing (not equal). 

M » -1 or -2. First or second order if VARI is strictly decreasing (not equal). 

N The number of points in the table (an integer), 

VARI The name of a one-dimensional array which contains the N values of the independent variable. 

VARD The name of a one-dimensional array which contains the N values of the dependent variable. 

Restrictions : All the numbers must be floating point. The values of the independent variable x in the 
table must be strictly increasing or strictly decreasing. The following arrays must be dimensioned by 
the calling program as indicated: VARI(N), VARD(N). 

Accuracy : A function of the order of interpolation used. 

References : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. The Macmillan Co,, c.1956, pp. 87-91. 

(b) Milne, Mlliam Edmund: Numerical Calculus. Princeton Univ. Press, c.1949, pp, 69-73. 

Storage; 430g locations. 

Error condition ; If the VARI values are not in order, the subroutine will print TABLE BELOW OUT OF 
ORDER FOR FTLUP AT POSITION xxx TABLE IS STORED IN LOCATION xxxxxx (absolute). It then 
prints the contents of VARI and VARD, and STOPS the program. 

Subroutine date ; September 12, 1969. 


no 
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